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IMAGE  SAMPLING  AND  INTERPOLATION 


I.  INTRODUCTION 

This  report  records  the  results  of  a  3-month  study  during  the  summer  of  1980  made 
by  the  author  while  on  temporary  assignment  to  the  Office  of  Naval  Research.  The  study 
was  begun  at  the  suggestion  of  Dr.  William  J.  Condell,  Jr.,  Head  of  the  Physics  Program  at 
ONR. 

This  study  includes  a  review  of  spline  functions,  in  particular  B-splines,  and  a  consid¬ 
eration  of  their  use  for  the  foveal  sampling  of  images.  Such  a  study  is  of  interest  to  the 
Navy.  In  many  applications  of  image  sampling,  as  for  example  in  a  missile  warhead,  only  a 
very  few  samples  can  be  handled.  Thus  methods  of  taking  these  samples  in  a  more  efficient 
manner  are  necessary.  It  was  thought  that  B-splines  might  be  useful  because  they  may  be 
used  for  interpolating  data  not  sampled  at  uniformly  spaced  points.  This  study  confirms 
this  idea  but  indicates  that  even  when  using  splines,  uniformly  spaced  data  make  numerical 
interpolation  more  efficient. 

Because  of  the  limited  time  available  and  the  newness  of  the  material  on  B-splines  to 
the  author,  all  of  the  possibilities  could  not  be  fully  developed.  In  particular  no  actual  nu¬ 
merical  work  with  images  was  done.  Without  actual  experience  with  images  no  study  of 
image  sampling  and  interpolation  can  be  taken  very  far.  The  scope  of  this  report  is  therefore 
limited.  In  it  the  author  reviews  the  well-known  Nyquist  sampling  theorem  and  points  out 
its  drawbacks  in  Section  II  in  order  to  motivate  the  use  of  splines.  In  Section  III  the  basic 
properties  of  B-splines  are  introduced  in  a  very  different  manner  than  is  usual  in  the  litera¬ 
ture.  While  this  material  lacks  rigor  it  should  be  much  simpler  for  the  nonmathematician  to 
follow.  The  use  of  the  so  called  pp-representation  and  B-representation  are  presented  in 
Section  IV  as  extensions  of  the  familiar  cubic  spline  methods.  And  finally  in  Section  V  a 
very  brief  account  is  given  as  to  how  all  of  this  might  be  applied  to  image  sampling  and  the 
foveal  problem.  The  conclusions  are  briefly  stated  in  Section  VI. 


II.  SAMPLING  AND  INTERPOLATION  OF  BANDLIMITED  IMAGES 

An  image  is  defined  in  this  report  as  a  two-dimensional  distribution  of  intensity  de¬ 
scribed  by  the  function  I(xy).  This  intensity  distribution  may  be  a  primary  image  obtained 
from  light  scattered  directly  from  a  scene  say  as  in  the  back  focal  plane  of  a  camera,  or  it 
may  be  a  secondary  image  formed  by  light  using  an  electronic  or  photographic  recording 
of  the  primary  image. 

If  the  function  I{xy)  has  finite  support  such  that 

I(xy)  *  0  if  |x|  >  a  or  |y|  >  b  ,  (2.1) 
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then  the  image  is  bounded.  All  images  formed  by  practical  instruments  are  bounded.  If, 
on  the  other  hand,  the  function  I(xy)  has  a  Fourier  transform 


n£,Tj)  = 


JJl(xy)e  2’n‘($x+riy  )dxdy  , 


(2.2) 


which  has  finite  support,  i.e., 


/'({,!?)  =  0  if  |£|>  a  or  |r?|>/3,  (2.3) 

then  the  image  is  bandlimited.  By  the  Paley-Wiener  theorem  of  Fourier  analysis  an  image 
cannot  be  both  bounded  and  bandlimited  [1] .  However,  this  theorem  reflects  a  basic 
difficulty  in  the  simple  mathematical  model  (or  in  achieving  ideal  optical  instruments  — 
depending  on  your  point  of  view)  because  images  do  appear  to  be  both  bounded  and  band- 
limited. 

Images  are  bounded  because  all  imaging  systems  contain  apertures  which  block  light 
outside  of  some  finite  area  in  the  image  plane.  Images  are  bandlimited  because  all  imaging 
systems  have  a  finite  resolving  power  and  can  in  no  case  resolve  image  detail  smaller  than 
about  a  wavelength  of  light.  In  the  following  we  assume  that  an  image  is  either  bandlimited 
or  that  it  is  bounded  but  not  both  so  that  we  may  use  this  mathematical  model. 

We  will  treat  the  case  of  a  linear  stationary  imaging  system.  Thus  a  bounded  image  can 
be  represented  by 

oo 

/(x,y)  =  f. f  h(x‘  -  *,/-  y)  0  (jc  o'')  dx'dy'  if  |x|  <a  and  |y|  <  b,  (2.4) 


=  0  otherwise, 


where  h(xy)  is  the  spread  function  (the  image  of  a  point),  and  0(x,y)  is  an  unbounded  and 
unbandlimited  “perfect”  image.  In  Eq.  (2.4)  I(x ,y)  is  bounded  but  not  bandlimited.  We  can 
represent  a  bandlimited  image  by 


OO 


h(x'  -  x,y'-  y)  O'(x\y')dx'dy'for  all  x  and  y  , 


(2.5) 


where 


0'(x',y)  =  0(x',y)  if  |x|  <  a  and  |y|  <  b  , 

=  0  otherwise. 

Now  I\xy)  represents  an  image  that  is  bandlimited  but  not  bounded.  In  the  rest  of  this 
section  we  will  consider  the  bandlimited  image  described  by  I'(xy). 
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To  record  an  image  it  is  often  convenient  (or  necessary)  to  sample  it.  For  example, 
a  vidicon  samples  an  image  continuously  along  horizontal  raster  lines  so  that  a  two-dimen¬ 
sional  spatial  function  I(xy)  can  be  represented  by  a  one-dimensional  temporal  function 
l '(f),  the  video  signal.  As  a  second  example,  a  photographic  emulsion  contains  silver  halide 
crystals  which  collect  light  and  are  converted  by  development  into  free  silver.  Thus  the 
crystals  sample  the  image  intensity  I(xy)  in  a  complicated  manner.  An  image  can  be  sam¬ 
pled  in  many  other  ways  that  can  be  selected  at  will  by  the  instrument  designer.  A  con¬ 
venient  choice  for  many  applications  [2]  is  to  sample  the  image  over  a  rectangular  mesh,  i.e., 

_  N.M 

£  I\xj)8(x- ndx)8(y- mdy)  (2.6) 

n=- N 
m=- M 


so  that  the  sampled  image  is  represented  by  the  discrete  data  values  I  where 

1  (xy )  =  Inm  if  x  =  ndx,y  =  mdy  ,  (2.7) 

=  0  otherwise, 


wh»  ic  dx,  dy  are  the  constant  sampling  intervals,  and  where  6  represents  the  Dirac  delta 
function.  Data  sampled  in  this  way  are  often  referred  to  as  gridded  data  in  the  interpola¬ 
tion  literature,  and  we  will  use  this  expression  in  this  report. 

To  record  or  transmit  the  sampled  image  we  need  only  deal  with  the  sequence  of  real 
numbers  Inm  rather  than  the  continuous  function  I'(xy y).  This  data  reduction  is  essential 
for  many  operations  on  the  image  as,  for  example,  digital  processing  by  a  computer. 

It  is  critically  important  to  the  accurate  and  efficient  representation  of  the  image  that 
the  sampling  intervals  have  the  values 


1 

d  =  —  and  d  =  —  .  (2.8) 

x  2<x  y  2(3  ’ 


These  equations  are  equivalent  to  the  Nyquist  criterion  which  states  that  a  bandlimited 
signal  must  be  sampled  at  uniformly  spaced  points  separated  by  half  the  period  of  the 
highest  Fourier  component  in  the  spectrum  of  the  signal.  It  is  shown  in  Appendix  A  that 
if  the  values  of  dx  and  dy  are  larger  than  specified  by  Eq.  (2.8),  then  aliasing  errors  occur 
which  prevent  accurate  recovery  of  the  image  from  the  sampled  data.  If,  on  the  other  hand, 
dx  and  dy  are  made  smaller  than  specified  by  Eq.  (2.8),  then  the  extra  samples  obtained 
are  wasted  because  no  additional  information  about  the  image  is  obtained. 


If  the  image  is  bandlimited  and  properly  sampled,  then  as  shown  in  Appendix  A  the 
image  can  be  recovered  exactly  by  using  the  formula 


N.M 

m-  M 


sin  d*(*-n<gj  sin  ~(y  rrdy) 


(x  -  ndx ) 


(y  -  mdy) 


(2.9) 
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We  will  call  this  sampling  and  interpolation  procedure  Nyquist  sampling.  Nyquist  sampling 
allows  the  image  I  (x,y)  to  be  recovered  at  every  point  from  just  the  sampled  values  Inm. 
Clearly  some  arbitrary  function  f(x,y)  contains  much  more  information  than  any  discrete, 
finite  set  of  sampled  values  fnm .  Therefore  we  must  conclude  from  Eq.  (2.9)  that  the  con¬ 
dition  that  /  (*,y)  be  bandlimited  is  a  very  strong  condition  which  greatly  limits  the  infor¬ 
mation  content  of  this  function. 

Often  it  is  not  possible  to  properly  sample  an  image.  Usually  the  number  of  samples 
that  car.  be  *  rken,  stored,  and  perhaps  transmitted  in  some  fashion  is  very  limited.  So  if  the 
image  is  not  properly  bandlimited  before  sampling,  aliasing  problems  develop.  This  well- 
known  problem  is  discussed  in  Appendix  A  and  has  been  examined  experimentally  [3] .  It 
appears  that  a  badly  aliased  image  is  usually  not  useful.  Therefore  a  completely  different 
approach  to  image  sampling  and  interpolation  is  needed  if  the  original  image  cannot  be 
properly  bandlimited.  Thus  the  requirement  that  an  image  be  properly  bandlimited  is  a 
serious  limitation  for  Nyquist  sampling. 

A  second  limitation  to  this  sampling  and  interpolation  method  is  that  it  is  restricted 
to  gridded  data.  It  is  often  more  useful  to  employ  a  limited  number  of  possible  samples  in 
a  more  efficient  manner  as  is  done  by  the  human  eye.  In  the  eye  samples  are  taken  much 
more  closely  together  over  a  small  region  near  the  center  of  the  field-of-view  called  the 
fovea  than  they  are  elsewhere.  The  result  is  that  we  have  a  much  better  ability  to  resolve 
fine  detail  in  the  foveal  image  than  elsewhere  in  our  field  of  view.  Samples  are  taken  so  far 
apart  near  the  edges  of  the  field  of  view  that  we  can  detect  little  more  than  motion  there. 

In  this  manner  our  eyes  make  the  best  use  of  a  limited  number  of  photoreceptors. 

In  this  report  we  are  particularly  concerned  with  optimizing  the  use  of  a  limited  num¬ 
ber  of  samples.  To  do  this  using  Nyquist  sampling  it  is  necessary  to  divide  the  image  into 
square  subareas  and  to  take  gridded  data  from  each  area.  The  sampling  intervals  may  be 
allowed  to  differ  from  one  area  to  another  but  must  be  the  same  within  each  area.  However 
there  remains  a  serious  problem.  It  is  very  difficult  to  individually  bandlimit  each  subarea 
of  the  image  so  that  the  Nyquist  criterion  is  satisfied  there.  Usually  the  entire  image  must  be 
bandlimited  in  the  same  way.  Then  if  the  image  is  properly  bandlimited  for  the  most  finely 
sampled  area,  it  will  be  aliased  elsewhere,  or  if  the  image  is  properly  bandlimited  for  the 
least  finely  sampled  area,  the  extra  samples  elsewhere  will  be  wasted  and  the  interpolated 
image  would  lie  uniformly  poor.  Thus  Nyquist  sampling  is  not  usually  of  much  use  under 
these  conditions. 

If  we  cannot  bandlimit  an  image  properly  so  that,  over  every  region,  the  image  satisfies 
the  Nyquist  criterion  locally,  then  Nyquist  sampling  should  not  be  used.  Fortunately  there  is 
a  useful  alternative.  It  is  always  possible  to  sample  an  image  at  any  set  of  points  and  then  to 
reconstruct  an  estimate  of  the  image  from  these  data.  It  is  not  possible  in  general  to  recon¬ 
struct  the  original  image  precisely  as  it  is  by  the  use  of  Nyquist  sampling;  however,  it  is  often 
possible  to  construct  a  useful  approximation. 

A  common  method  for  reconstructing  a  one-dimensional  function  f(x)  from  scattered 
data  (samples  taken  at  random)  is  to  set 

fix)  =  fn  ,  where  x  =xn  ,  (2.10) 
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so  that  at  the  original  sample  points  xn  the  function  f(x)  has  its  original  values  and  then  to 
use  a  draftsman’s  spline  (a  tool  like  a  rubber  ruler)  to  continue  the  data  smoothly  in  be¬ 
tween.  This  gives  a  much  better  representation  of  the  original  function  than  does  Nyquist 
sampling  with  serious  aliasing.  As  a  practical  matter  for  a  two-dimensional  function  it  is 
better  to  employ  the  numerical  methods  of  fitting  surfaces  to  scattered  data.  Here  we  rep¬ 
resent  the  draftsman’s  tool  by  a  piecewise  polynomial  “spline”  function  which  is  used  to 
connect  the  data  points  in  a  satisfactory  manner. 

In  the  following  sections  we  introduce  some  of  the  fundamentals  of  spline  interpola¬ 
tion  and  discuss  the  application  of  these  methods  to  image  interpolation. 


III.  B-SPLINE  FUNCTIONS 

fl-splines  or  basis  splines  were  first  introduced  by  Schoenberg  [4]  as  a  set  of  functions 
defined  by 


. . S  £(*?)•->"> 


(3.1) 


For  k  =  1  it  is  well  known  to  mathematical  physicists  that 

^1(0  =  F7’(^^-))  =  Rect(f) 


=  i  if  in  <  - 
2 

=  0  otherwise, 

where  FT  represents  a  Fourier  transform.  It  is  also  well  known  that  by  use  of  the  con¬ 
volution  theorem  Eq.  (3.1)  can  be  rewritten 


(3.2) 


FT 

~  /  sin(£/2)  \ft-f  ' 

®FT 

l sin(£/2)  V 

_\  €/2  ) 

.  V  m  ) . 

(3.3) 


so  that  Mk(t)  can  be  built  up  from  the  Rect  ( t )  function  via  k  -  1  convolutions.  It  is  clear 
from  this  that  because  Rect  {x)  has  small  support  (i.e.,  vanishes  outside  of  the  domain 
- 1/2  =  t  5?  1/2)  that  all  of  the  Mk(t)  built  up  by  repeated  convolution  with  Rect  (f)  will 
similarly  have  small  support  (i.e.,  vanish  outside  of  the  domain  -k/2  ^  t  S  k/2).  This  is  one 
of  the  most  useful  properties  of  the  B-splines  ([5] ,  p.  109). 

The  B-splines  up  to  cubic  order  are  given  in  Fig.  3.1  for  the  case  of  uniform  sampling. 
Note  that  Mk(t)  is  a  peacewise  polynomial  (pp )  function  of  order  k  -  1.  It  is  easily  demon¬ 
strated  that  is  continuous  with  continuous  derivatives  up  to  order  k  -  2.  Higher  order 

derivatives  are  continuous  everywhere  but  at  the  break  points  between  polynomial  seg¬ 
ments.  These  break  points  are  called  knots  in  the  literature  on  B-splines. 
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The  definition  of  B-splines  given  by  Eq.  (3.1)  holds  only  for  the  special  case  of  an 
equispaced  knot  sequence,  that  is  for  polynomial  segments  of  equal  length  as  shown  in 
Fig.  3.1.  However,  Curry  and  Schoenberg  [6]  showed  that  if  the  definition  of  the  B-splines 
is  expressed  in  terms  of  the  feth  order  divided  difference  as  described  here  in  Appendix  B, 
then  the  definition  can  be  extended  to  an  arbitrary  knot  sequence.  This  admits  a  much  more 
general  class  of  pp-spline  functions  composed  of  polynomial  segments  of  arbitrary  length 
and  with  discontinuities  of  arbitrary  order  at  the  break  points  between  segments. 

To  extend  B-splines  to  allow  arbitrary  knot  sequences  [fj,  t2  ...  4]  we  redefine  the 
B-spline  of  order  1  at  location  i  by 

Bj.iU)-1  if  , 


=  0  otherwise,  (3.4) 

which  agrees  with  the  definition  in  Eq.  (3.1)  and  then  generate  higher  order  B-splines  using 
the  recursion  relation  (see  Appendix  B,  Eq.  (B.17)) 


— T lB^- 


i+k-  1 


i<')+r 

(+ 


w  -  * 


k  ~  ii+i 


B; 


i+l,fc 


-!(*>• 


(3.5) 


This  is  equivalent  to  the  definition  of  B,  k(t)  in  the  mathematical  literature  which  depends 
on  the  rather  complicated  concept  of  the  divided  difference  (see  Appendix  B).  If  the  knot 
sequence  ( f x ,  <2  ...  fg]  is  monotonically  increasing  with  constant  interval  between  knots, 

Eq.  (3.5)  is  equivalent  to  the  definition  in  Eq.  (3.1).  Otherwise  it  is  not. 

B-splines  are  very  important  in  the  theory  of  pp-functions  because  they  form  a  basis 
for  such  functions  (hence  the  name  basis  splines).  Any  pp-function  can  then  be  represented 
as  a  superposition  of  B-splines  of  the  form  ([5] ,  pp.  119-120): 

n 

4(i)=E«A,k(‘)  0.5) 

;=  i 

or  if  tj  <  t  <  £y+  j  by 


4(0=  £  «B(,fe(0  (3.6) 

i  =  j-  fe  +  1 

where  4(0  is  a  PP  spline  function  of  order  k  for  the  arbitrary  knot  sequence  [4,  £2  •••  £g] , 
and 


n  =  fe«- 


i=2 


(3.7) 
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The  knots  are  a  sequence  of  points  between  the  polynomial  segments  of  the  B-splines. 

They  may  be  coincident  as  indicated  by  the  parameter  tq(<fe)  in  such  a  manner  that  at 
some  point  t  there  are  ( k  -  i>-)  knots  t,.  To  use  a  pp-function  to  interpolate  sampled  data 
the  knot  may  represent  the  location  of  data  samples.  They  may  be  separated  by  arbitrary 
intervals  corresponding  to  nonuniform  sampling  intervals.  The  coincident  knots  are  useful 
for  specifying  the  number  of  continuity  conditions  the  pp-function  satisfies  over  the  knot 
point.  At  a  point  with  a  single  knot  the  pp-function  and  its  derivatives  up  to  the  order  k  -  1 
must  be  continuous  over  the  point,  but  at  a  point  with  k  -  knots  the  pp-function  will 
satisfy  only  v •  continuity  conditions  over  the  point.  Multiple  knots  are  also  useful  for 
osculatory  interpolation  in  which  we  specify  data  concerning  derivatives  as  well  as  values  of 
the  sampled  function. 

Interpolation  and  the  use  of  B-splines  are  discussed  in  the  following  sections. 


IV.  INTERPOLATION  IN  ONE  DIMENSION 

The  interpolation  of  data  sampled  along  a  single  line  is  introduced  in  this  section.  The 
more  complicated  two-dimensional  interpolation  problem  which  applies  to  images  is  dis¬ 
cussed  in  the  next  section. 

The  hest  known  and  probably  the  most  commonly  used  spline  interpolation  method 
employs  data  fitting  with  piecewise  polynomial  cubic  splines.  The  application  of  this 
method  in  one  dimension  is  clearly  described  in  many  standard  works  on  numerical  analysis 
([7] ,  p.  474-491).  We  will  describe  this  method  first  as  an  introduction  to  the  more  general 
methods. 

For  the  usual  pp-cubic  spline  interpolation  a  different  cubic  polynomial,  for  example 

f(x)  =  a,(.r  -  x,)3  +  6(x  -  x;)2  +  e((x  -  *,)  +  d  .  (4.1) 

is  fitted  between  each  point  x((  1  <  i  <■.  n )  at  which  data  have  been  sampled  and  the  next 
larger  point  xj+  j .  The  sample  points  become  break  points  between  these  polynomial  seg¬ 
ments.  The  coefficients  in  each  segment  are  determined  by  making  restrictions  on  the 
pp-cubic  spline.  The  most  usual  conditions  are  that  the  pp-function: 

(1 )  is  continuous  over  the  break  points  (n  -  2  equations) 

(2)  gives  the  sampled  data  value  at  each  of  the  break  points  ( n  equations) 

('3)  has  continuous  first  and  second  derivatives  over  the  break  points  (2(n  -  2) 
equations). 

Thus  for  the  n  sampled  data  values  we  obtain  4n  -  6  equations  for  the  4  ft  -  4  coefficients 
a,,  br  rr  dt  in  Kq.  (4.1).  Clearly  we  require  two  additional  conditions  in  order  to  specify  the 
coefficients.  These  can  be  made  arbitrarily,  for  example,  by  specifying  the  derivatives  of  the 
pp-function  of  each  of  the  two  end  points.  The  An  -  4  equations  are  easily  reduced  to  n 
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equations  analytically,  and  solution  of  the  remaining  n  equations  containing  the  sampled 
data  values  reduces  to  inverting  an  n  X  n  tridiagonal  matrix  using  a  computer. 

There  are  several  variations  to  pp-cubic  spline  interpolation  that  can  sometimes  give  a 
more  satisfactory  result  in  some  special  cases.  One,  of  course,  is  to  choose  differently  the 
two  arbitrary  restrictions  described  in  the  last  paragraph.  Another  is  to  relax  the  conditions 
that  the  second  derivative  be  continuous  over  the  break  points  and  instead  require  that  the 
first  derivative  be  specified  at  the  break  points  by  sampled  data.  Clearly  there  is  a  lot  of  flex¬ 
ibility  to  pp-interpolation.  This  is  brought  out  more  clearly  in  the  generalized  methods 
which  we  will  now  briefly  introduce. 

The  pp-cubic  spline  interpolation  can  be  easily  modified  by  using  polynomials  of 
higher  or  lower  order.  Higher  order  polynomials  of  order  k  (cubic  is  k  =  4)  can  be  used  to 
give  a  pp-interpolation  function  which  has  continuous  derivatives  to  higher  order  k  -  2  over 
the  break  points  but  at  the  cost  of  more  computation.  The  matrix  is  still  n  X  n  but  is  no 
longer  tridiagonal  but  has  k  nonzero  semidiagonals  instead  (fe-diagonal).  Lower  order  poly¬ 
nomials  greatly  reduce  the  computations  but  give  less  smooth  pp-functions. 

Once  the  interpolation  calculations  have  been  carried  out,  the  pp-function  is  stored  in 
ihe  computer  for  use  in  calculating  values  of  the  function  f(x)  which  was  represented  orig¬ 
inally  by  the  sampled  data.  To  store  a  pp-representative  for  some  function  f(x)  in  a  computer 
we  require  t  (5] ,  p.  88) 

(1)  the  integers  k  and  n  giving  the  order  and  number  of  the  polyn  omial  pieces, 
respectively. 

(2)  the  strictly  increasing  sequence  x} ,  x2,  x3,  ...  x^  +  J  of  its  breakpoints  (which  do 
not  in  general  have  to  be  equally  spaced),  and 

(3)  the  matrix 


C»  - 


y 


r  1 


dxi- 


f{x)  for  1  </  <  k,  and  !<('<« 


of  its  (right)  derivatives  at  the  breakpoints.  Then  the  pp-function  can  be  used  to  represent 
the  sampled  function  over  each  domain  x,  =  x  =  xl+  j  by 


flx)  = 


k-i  (x -$,.)"" 

S  Cm  + 1  ,i  m  i 
m  =0 


(4-2) 


This  is  nothing  but  a  generalization  from  the  cubic  spline  representation  to  splines  of  k  order 
where  for  k  =  4  we  have  a,  =  C4|/ 6,  6,  =  C3j/2,  c,  =  C2i,  d,-  =  C1(  so  that  Eqs.  (4.1)  and  (4.2) 
become  the  same.  Clearly  the  use  of  higher  order  splines  provides  greater  flexibility  for 
improving  data  interpolation  but  at  the  cost  of  greater  computation. 
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Still  greater  flexibility  can  be  obtained  by  the  use  of  an  entirely  different  representa¬ 
tion  based  on  the  B-splines  introduced  in  the  last  section.  In  interpolating  a  sequence  of 
data  points,  B-splines  (for  example  as  show  a  in  Fig.  3.1,  if  the  data  points  are  equally  spaced) 
of  some  fixed  order  k  are  fitted  over  k  successive  knots.  Each  of  the  knots  (except  the  free 
knots  discussed  later)  can  represent  a  sampled  data  point.  The  B-splines  unlike  the  poly¬ 
nomial  segments  described  above  are  allowed  to  overlap  in  such  a  manner  that  a  new  B-spline 
begins  at  each  knot  location.  Thus  we  have  a  rather  different  sort  of  interpolation  method 
with  additional  redundancy  which  leads  to  greater  flexibility  in  application. 

The  calculation  of  a  B-representation  for  a  sampled  function  f(x)  reduces  to  the  inver¬ 
sion  of  a  matrix  by  a  computer.  The  matrix  is  again  fc-diagonal;  however,  it  is  now  of  much 
higher  order  due  to  the  extra  coefficients  required  to  specify  the  over'apping  B-splines.  The 
order  of  the  matrix  is  now  n  X  n  where 

n  =  fe«  -  E  ",  .  (4.3) 

i 

where  k  is  the  order  of  the  B-splines,  t  is  the  number  of  actual  data  samples,  and  Vj  is  the 
order  of  continuity  over  the  knot  points  (i.e.,  vi  =  0,  no  continuity;  i>  =  1,  continuous  f(x) 
only;  =  2,  continuous  f(x)  and  first  derivative  of  f(x);  etc.). 

Once  the  B-representation  has  been  calculated  from  the  sampled  data  it  is  stored  in  the 
computer.  Then  the  B-representation  for  some  originally  sampled  function  f(x)  consists  of 
([51 ,  p.  119): 

ill  the  integers  k  and  n  (n  =  kV  -  i>()  giving  the  order  of  the  B-spline  segments  in  the 
interpolating  function  and  the  number  of  linear  parameters  required  to  represent  the 
function, 

(2)  the  vector 


<1  <  *2  <  f3  <  •  <  **  <  *1  <  <  tn  +  X  (4'4> 

containing  the  knots  (possibly  partially  coincident)  in  increasing  order,  and, 

(3)  the  vector  a-,  1  <  /  <  n  of  the  coefficients  of  f  with  respect  to  the  B-spline  basis 
as  shown  in  Eq.  (3.6). 

The  knots  to  +  1  can  be  actual  sample  points  whereas  the  free  knots  tx  to  tb  and  <„  +  1 
to  tn+l,  are  arbitrary  and  can  be  chosen  to  optimize  the  interpolation.  The  proper  choice  of 
knots  is  a  rather  complicated  subject  ([5] ,  Chap.  XIII)  that  we  will  not  attempt  to  cover  at 
the  level  of  this  report.  A  proper  understanding  of  this  subject  gives  the  computer  program¬ 
mer  great  flexibility  in  controlling  the  properties  of  the  interpolating  function. 
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Once  the  calculations  have  been  carried  out  and  the  B-representation  data  stored  in  the 
computer  then  the  value  of  the  sampled  function  f(x)  can  be  calculated  between  each  pair 
of  knots  Xj  <x  <x;+1  by  use  of  Eq.  (3.6),  i.e.. 


f(X)  =  £  a,B,  k(x)  . 
i=j~  k  +  1 


(4.5) 


In  this  manner  we  can  calculate  a  value  of  f(x)  for  each  point  within  the  domain  bounded 
by  the  data  samples.  At  the  sample  points  we  recover  the  sampled  values  and  in  between  we 
get  a  smoothly  continued  set  of  values  which  are  determined  by  the  data  and  by  the  knots 
chosen  by  the  programmer. 

If  the  original  sampled  function  contains  noise  such  that 

f(x)  =  g(x)  +  n(x)  (4.6) 

where  g{x)  is  the  signal  and  n(x)  is  the  noise,  then  spline  interpolation  of  samples  from  f(x) 
can  often  be  chosen  such  that  the  resulting  estimate  of  f(x)  is  smoothed.  For  this  we  choose 
the  interpolation  knots  to  be  at  the  sample  points,  and  we  employ  spline  interpolation  with 
an  auxiliary  condition  that  the  pp-function  minimizes  a  roughness  expression. 

In  cubic  spline  interpolation,  for  example,  to  smooth  the  sampled  data  we  replace  the 
usual  conditions  with  the  condition  that  for  n  sample  points  the  pp-function: 


(1)  is  continuous  over  the  break  points  (n  -  2  equations), 

(2)  gives  values  a,  at  the  break  points  (n  equations), 

(3)  gives  continuous  first  derivatives  at  the  break  points  (n  -  2  equations), 

(4)  gives  zero  second  derivatives  at  the  end  points  (2  equations), 

(5)  gives  values  c,-  for  the  second  derivatives  at  the  break  points  (n  -  2  equations). 

This  gives  the  system  of  4n  -  4  equations  required  to  determine  the  coefficient  in  Eq.  (4.1). 
The  2 n  parameters  a,  and  c,  are  determined  by  requiring  that  the  roughness  expression 
((5J ,  p.  238) 


Rj(p)  =  p  Y, 
1=  1 


5y, 


+ 


4(1 -p) 


n~  1 

£  **i(  C?+C«C.+  I  +chl  )  (4.7) 

i=i  ' 


be  minimized  for  some  value  of  p  (0  <  p  <  1),  where  y,  is  the  sampled  value  at  the  *,■  point, 
6y(  is  an  estimate  of  the  variance  of  y,-,  and  Ax,  =  xi+  j  —  jcf.  If  we  set  p  =  0,  the  pp-function 
is  optimally  smooth  but  inaccurate,  but  if  we  setp  =  1,  the  pp-function  is  accurate  (i.e., 
a,  =  y,)  but  not  smooth.  In  practice  a  value  of  p  is  chosen  to  give  the  best  trade-off  in  the 
opinion  of  the  programmer. 
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In  the  next  section  we  consider  the  extension  of  one-dimensional  interpolation  to 
images. 


V.  IMAGE  INTERPOLATION 

The  extension  of  the  interpolation  methods  introduced  in  the  last,  section  from  one  to 
two  dimensions  is  not  trivial  except  for  gridded  data.  We  are  particularly  interested  in  inter¬ 
polating  scattered  data  in  order  to  implement  foveal  sampling.  That  is,  we  would  like  to 
interpolate  data  which  were  sampled  at  points  in  the  original  image  that  were  much  closer 
together  over  the  foveal  region  than  outside  of  it  as  discussed  in  Section  II.  Thus  we  will 
consider  scattered  data  first. 

Consider  an  image  which  has  been  sampled  in  the  following  manner.  Let  a  general 
curvilinear,  orthogonal,  two-dimensional  coordinate  system  span  the  image  plane  which  is 
assumed  to  be  finite  but  of  arbitrary  shape.  The  coordinate  curves  must  describe  a  pattern 
of  lines  which  intersect  at  right  angles  as  shown  in  Fig.  5.1.  The  data  samples  are  taken  from 
the  image  at  points  £  =  £,  along  curves  of  constant  17  =  17,.  We  take  £,,  r?,  to  be  closely  spaced 
in  the  foveal  region  and  more  widely  spaced  elsewhere.  An  example  of  a  set  of  such  sam¬ 
pling  points  is  shown  in  Fig.  5.2.  In  this  example  only  80  samples  are  shown;  however,  to 
sample  this  area  with  gridded  data  and  at  the  foveal  sampling  rate  would  require  144 
samples.  We  have  reduced  the  number  of  samples  44%  and  maintained  the  foveal  resolution, 
but  at  a  cost.  We  cannot  precisely  recover  the  image. 

To  recover  an  approximation  to  the  image  by  using  the  sampled  data  we  employ  one¬ 
dimensional  spline  interpolation  first  along  the  curves  of  constant  £  =  £,  and  then,  using  the 
interpolated  values,  along  all  required  curves  of  constant  ij.  In  this  manner  an  estimate  of 
/(£,  i?)  at  any  point  over  the  image  can  be  recovered.  Since  the  accuracy  of  the  estimate  at 
a  point  will  improve  as  the  point  nears  a  sampling  point,  the  approximation  to  the  image 
will  be  more  accurate  over  the  foveal  region  where  the  sampling  points  are  closer  together. 
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Fig.  5.2  —  Example  of  foveal  sampling 


Many  methods  have  been  invented  for  interpolation  of  scattered  data  in  two  dimen¬ 
sions.  A  general  review  article  on  this  subject  was  published  by  Schumaker  [8] .  The  general 
problem  of  two-dimensional  spline  interpolation  has  not  been  developed  as  well  as  in  one 
dimension.  It  is  always  possible  to  use  a  two-dimensional  pp-function  over  each  polygon¬ 
shaped  region  bounded  by  nearest  sample  points.  In  general,  the  calculation  of  coefficients 
representing  the  pp-function  cannot  be  computed  very  efficiently.  The  methods  used  de¬ 
pend  on  the  geometry  of  the  sampling  points.  Kor  gridded  data  the  pp-function  reduces  to 
the  tensor  product  of  two  one-dimensional  pp-functions  and  the  calculation  of  the  coef¬ 
ficients  becomes  particularly  efficient.  Thus  there  is  a  distinct  advantage  to  be  had  using 
gridded  data  ([5] ,  p.  332,  [8] ,  p.  218). 

To  take  advantage  of  the  numerical  efficiency  of  spline  interpolation  of  gridded  data 
and  at  the  same  time  make  more  efficient  use  of  a  limited  number  of  samples  we  will 
consider  again  the  foveal  sampling  approach  described  in  Section  II.  Let  the  image  plane  be 
decomposed  into  rectangular  subareas.  Each  area  is  sampled  using  gridded  data,  but  the 
sampling  intervals  dx  and  dy  in  Eq.  (2.6)  can  be  chosen  differently  in  different  subareas. 
Over  the  subareas  where  high  resolution  will  be  required  we  make  dx  and  dy  small,  whereas 
over  regions  where  lower  resolution  is  adequate  we  make  dx  and  dy  proportionally  larger. 
Since  we  intend  using  spline  interpolation,  the  image  need  not  be  properly  bandlimited 
within  each  subarea. 

An  estimate  of  the  image  is  recovered  from  the  samples  in  each  subarea  by  the  use  of 
surface  spline  functions  given  by 

Bfjixy)  =  Bi  k(x)Bj  k(y)  ,  (5.1) 
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where  B[  ^(.v)  and  ;.(y)  are  defined  in  Eq.  (3.5)  ([8] .  Eq.  (3.30)).  These  splines  are  super¬ 
imposed  to  yield  a  pp-function  interpolating  the  data  as  given  by 

/(-v,y)  =  Y,  <5~» 

i.j 


in  which  the  atJ  coefficients  are  to  be  computed  in  a  manner  very  similar  to  computation  in 
the  one  dimensional  case  of  the  last  section.  Computer  programs  exist  for  two-dimensional 
bicubic  interpolation  ([8] ,  p.  221).  The  computation  of  the  a,y  coefficients  can  be  carried 
out  more  efficiently  than  by  straightforward  extension  of  the  one-dimensional  method  to 
two  dimensions  ([5] ,  p.  343).  Some  simple  experiments  with  cubic  spline  interpolation  of 
image  data  were  done  by  Hou  and  Andrews  [9,10] ,  but  their  work  was  limited  to  equi- 
spaced  data  taken  at  the  knots.  Much  more  experimental  work  should  be  done  to  develop 
and  extend  this  as  a  useful  technique. 


VI.  CONCLUSIONS 

The  conclusions  reached  in  this  study  can  be  stated  very  simply.  The  best  method  of 
sampling  an  image  is  by  the  use  of  Nyquist  samples  provided  the  image  is  properly  band- 
limited  and  gridded  data  are  used.  If  the  number  of  samples  that  are  allowed  is  very  limited, 
then  the  best  method  is  to  divide  the  image  into  subareas  of  gridded  data.  If  the  image  can 
be  properly  bandlimited  within  each  subarea,  then  again  Nyquist  sampling  is  best,  but  if 
the  image  cannot  be  properly  bandlimited,  then  spline  interpolation  is  better. 

These  conclusions  do  not  hold  for  some  unusual  situations  like,  for  example,  the  ex¬ 
treme  case  of  absolute  minimum  samples  and  unlimited  computer  capability  where  the  more 
general  use  of  scattered  data  and  general  spline  interpolation  is  probably  nearer  to  optimum. 
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Appendix  A 

THE  SAMPLING  THEOREM 

Consider  a  continuous  function  fix )  defined  for  every  x  within  the  domain 
-oo  <  .v  We  can  represent  the  sampling  of  this  function  hy  multiplying  fix)  by  a  Dirac 

Comb,  i.e.. 


oO 

fix)  =  £  f(x)l)(x  -  nd) 

n  =  - 


( A.l ) 


where  d  is  the  distance  between  the  samples,  and  from  Eq.  (A.l)  we  see  that 

fix)  <*  f(x)  if  x  is  an  integer  multiple  of  d,  (A.2) 

=  0  otherwise. 


is  a  collection  of  samples.  This  operation  is  illustrated  by  Fig.  A.l. 


Kig  A.l  —  Sampled  data  fix) 

To  transmit  the  sampled  function  /  (x)  we  need  only  send  the  values  of  f  (x)  for 
x  =  nd  as  a  time-ordered  sequence  of  numbers.  At  the  receiving  end  we  want  to  reconstruct 
the  function  fix)  from  these  numbers  as  well  as  we  can.  Thus  we  require  the  inverse  of 
Eq.  (A.l)  which  gives  fix)  as  a  function  of  f  (x). 

To  invert  Eq.  (A.l)  we  first  take  its  Fourier  transform,  i.e., 

F(i>)  =  f  J2  fix)b{x  ~  nd)e2Tnvxdx 

J  -  oo  h  3-  -  oo 


oo  ^.co 

=  Yl  I  f{x)b{x  ^  nd)  e2nn’*dx  (A. 3) 

n  =  -  oo  •/-  oo 


oo 

=  Fiv)®  ^2  e2nivnd  ^ 

n  =  '  00 
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where  we  interchanged  the  orders  of  integration  and  summation,  then  used  convolution 
theorem  to  remove  /'(x)  from  the  kernel  of  the  integral  and  finally  carried  out  the  integra¬ 
tion  using  the  sifting  property  of  the  Dirac  delta  function.  Consider  the  summation  in 
Eq.  (A. 3)  as  a  function  of  t\  For  v  -  m/d  (rn  =  any  integer)  we  have 


QO 

£  e2nimn  _ oe 

n=-°° 


(A. 4) 


However  if  v  =£  m/d,  we  have  instead 


£ 


e2nivnd 


n  =  -  °° 


0, 


(A. 5) 


since  e2nix  averages  to  zero.  Thus,  we  have  from  Eq.  (A. 3) 


oo 

F(i’)  =  F(v)C-)  53  5(v-  m/d)  , 

m  -  -  00 


(A. 6) 


which  represents  the  frequency  spectrum  F(v)  of  the  original  function  f{x)  convolved  with 
another  Dirac  Comb  with  a  spacing  inversely  proportional  to  d.  This  is  shown  in  Fig.  A. 2.  If 
the  original  function  was  bandlimited  before  sampling  such  that  F(v)  =  0  outside  of  the 
domain  -  l/2d  <  l/2c/,  then  the  various  orders  in  this  figure  do  not  overlap  which  is  the 

case  as  shown.  Then  to  remove  the  effect  of  sampling  we  need  only  select  the  central  order 
by  multiplying  F(v)  by  a  Rect  function  in  the  manner 


F(v)  =  F(v)  Rect  (vd) 


=  F(v) 


(A. 7 ) 


=  0  otherwise  . 


Ei(!  A. 2  —  Spectrum  of  sampled  data  F(i’) 
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Equation  (A. 8)  shows  that  we  can  reconstruct  the  function  fix)  for  every  x  from  just  the 
sampled  values  at  x  =  nd  by  convolving  the  samples  with  a  sme  function  spline.  This  is  the 
usual  sampling  theorem.  If  the  original  function  fix  I  was  not  bandlimited  as  assumed  in 
Fig.  A. 2  so  that  the  orders  overlap  as  shown  in  Fig.  A. 3,  then  the  spectrum  of  fix)  is  not 
simply  repeated  as  before  but  high  frequencies  from  one  order  are  confused  as  lower  fre¬ 
quencies  by  the  next  order.  If  we  attempt  to  use  the  interpolation  formula  in  F'q.  (A. 8),  we 
obtain  the  function  f  ix )  with  the  spectrum  shown  in  Fig.  A. 4.  In  the  region  of  overlap  the 
frequencies  are  given  by 

F  O’)  =  Fiv)  +  Fit'  -  i\ld)\  .  (A. 9) 


CARTER 


The  contribution  from  the  second  term  in  Eq.  (A. 9)  represents  beats  between  the  sampling 
raster  and  the  frequencies  in  f(x).  This  appears  in  the  interpolated  function  f'(x)  as  errors 
termed  aliasing  errors.  The  phenomena  is  called  aliasing.  The  requirement  that  Fiv)  =  0 
unless  -  1/d  <  v  <  1/d  which  avoids  aliasing  is  equivalent  to  the  well-known  sampling  criteria 
that  f(x)  must  be  sampled  twice  per  period  of  the  highest  frequency  component  in  fix).  In 
two  dimensions  the  sampling  theorem  for  a  rectangular  array  of  equispaced  samples  can  be 
generalized  to  give 


sin 


f(xy)=  X)  f(nd,md)- 


r  7r  i 

r 

sin 

(y  -  mdy ) 


■J 


(x  - nd } 


(y  -  mdy) 


(A. 10) 


A  similar  procedure  has  been  worked  out  for  uniform  sampling  in  circular  coordinates  which 
is  rather  more  complicated  [A1,A2] . 
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Appendix  B 


B-SPLINES  REPRESENTED  AS  A  DIVIDED  DIFFERENCE 


There  is  another  representation  for  the  B-splines  which  is  often  more  useful  than 
Eq.  (3.1).  To  obtain  this  representation  consider  any  n  times  differentiable  function  h(x) 
and  consider  the  Taylor’s  series  representation 


h(x)  = 


n-  1 

£ 


i=0 


(X  -  g)1 
i 


(x  -  t)”'1 
(n-  1)! 


h("\t)dt  . 


(B.l) 


We  consider  a  sequence  of  points  on  x  (sometimes  called  a  knot  sequence)  given  by 
x, ,  x2,  ...  x„  in  which  x(  increases  monotonically  with  i  and  with  a  constant  interval 
Ax  =  x1+  j  -  Xj.  We  define  the  first  divided  difference  of  h(x)  at  x,  by 


[*/’X(+i  ]/»(')•■ 


h(xt)  -  h(xi+l) 
~~xi~xi* 1 


(B.2) 


and  similarly  the  kth  divided  difference  is  given  by  ([5] ,  p  8,  property  VIII) 

-  [•*>  xr-  lfXr+l’  ••  ~  [Xi>  •> 


[^,x1+1,...,xi+fe]h(.): 


x,  -  X . 


(B.3) 


where  xr  and  xs  are  any  two  different  points.  Thus  we  see  from  Eq.  (B.3)  that  a  feth  order 
divided  difference  of  h(x)  can  be  built  up  by  taking  k  first  order  divided  differences  of  h(x). 
The  divided  difference  is  clearly  very  closely  connected  to  the  derivative  of  the  function 
([5] ,  p.  8,  property  VII)).  If  we  take  the  nth  order  divided  difference  of  Eq.  (B.l),  we  find 
that 


(,)(-  -  t)"'1 

(»-  D! 


hln\t)dt 


where  we  have  used  the  property  ([5] ,  p.  6,  property  V) 

[*,,  ....  *;+„]  "  a)m  =0  if  m  <  n  , 


(B.4) 


(B.5) 
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which  can  be  proven  by  direct  application  of  Eq.  (B.3).  If  we  define 

(x-  t)""1  =  (x-  t)n~1,if(x-  t)>  0, 
=  0  otherwise, 

then  we  have  (for  b  >  x) 

(x  -  tf'1 


(B.6) 


Jr*  -  t)n~ 1  rb  (x  -  *)+  1 

1  7—  -h(n)(t)dt=  f  — — —  h("\t)dt 

a  (n~  1)!  Ja  (n  -  1)! 


(«-  1) 

and  by  substitution  into  Eq.  (B.4)  we  find  that 


(B.7) 


t*  -j u ^ — k(n)^dt- 


(B.8) 


Now  if  we  let 


h(x)  =  e'$* 

a  =  -  oo 
6  =  °°, 


Eq.  (B.8)  becomes 


(B.9) 


(B.10) 


[*,•  ■’Xi+n]e‘iX  =  <«'$>"  J ^  - ^7^ - e*  dt  • 

The  RHS  of  Eq.  (B.10)  is  just  the  Fourier  transform  of  the  nth  divided  difference  of 
(x  -  t)+~  1  and  the  LHS  can  be  evaluated  using  (B.3)  to  give 

[-. . 

when  Ax  =  xl  +  1  -  x,,  and  where  x  =  (xJ+„  +  x,  )/2  is  the  middle  point  of  the  sampling 
domain. 

By  substitution  from  (B.ll)  into  (B.10)  we  have 

"■"/> . *.■•.]  <B.H, 
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which  by  Fourier  inversion  becomes 


1  f  [  sin  ^/2)  " 
2rr  1  „  L  Z/2  . 


(B.13) 


where  we  have  rescaled  a:  so  that  Ax  =  1  and  *  =  0.  By  comparison  of  Eq.  (B.13)  with 
Eq.  (3.1 )  we  have 


M,M)  =  k[x Xi+k]  (•-  Of1 


(B.14) 


the  fcth  spline  is  the  A'th  divided  difference  of  the  function  l:(x  t)+  ^ .  I  his  is  the  relation 
that  is  almost  always  used  in  the  mathematical  literature  to  define  the  B-splines. 


Although  the  original  definition  of  /1-splines  as  given  in  Eq.  (3.1 )  holds  only  for  a 
sequence  of  equally  spaced  knots  the  definition  in  Eq.  (B.14)  can  be  generalized  to  an 
arbitrary  knot  sequence  ([5] ,  p.  108),  i.e.. 


*, •.*<*>  =  (w  -  r/)[rf  T/>  J 


(•  -  / 


>/.-  1 


(B.15) 


where  the  knot  sequence  r,,  r,+  j ,  are  arbitrary  points  on  a  line  (in  any  order,  with 

any  spacing,  and  possibly  with  repeated  values). 


By  applying  Leibniz’  formula  for  the  h th  divided  difference  of  a  product  to  the  partic¬ 
ular  product 

(t-x)*-'  =(/-*)«•■ xt2  (B-16) 


we  can  show  ([5] ,  p.  130)  that  the  generalized  B- splines  obey  the  recurrence  relation 


t  -  t. 


ti+k  -  1 


+  . — r~Bi* i 

1  V  li*k  li+ 1 


fr-  1 


(0 


(B.17) 


From  Eq.  (B.17)  and  the  evident  fact  (from  Eq.  (B.15))  that 

«,  ,(/)=  1  >f  <  t  <  t,-*  i  • 

=  0  otherwise,  (B.18) 

we  have  a  simple  way  of  generating  the  //-splines  in  a  computer.  Because  of  the  complexities 
of  dealing  with  divided  differences  we  use  Eqs.  (B.17)  and  (B.18)  to  define  generalized 
//  splines  in  this  report. 


